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ABSTRACT 

A new theory of gravitational shocking based on time-dependent perturbation 
theory shows that the changes in energy and angular momentum due to a slowly 
varying disturbance are not exponentially small for stellar dynamical systems in 
general. It predicts significant shock heating by slowly varying perturbations 
previously thought to be negligible according to the adiabatic criterion. The 
theory extends the scenarios traditionally computed only with the impulse 
approximation and is applicable to a wide class of disturbances. The approach 
is applied specifically to the problem of disk shocking of star clusters. 

1. Introduction 

A wide variety of astronomical problems require solving for the evolution of a bound 
stellar system in the external gravitational held of a larger embedding system. Examples 
include globular clusters on eccentric orbits in the Galaxy, infalling satellite galaxies, 
galactic star clusters in the disk, and the evolution of young associations in molecular 
clouds. 

In many cases, the characteristic time scale for the external perturbation is neither 
slower or faster than all orbital periods of the bound cluster, but somewhere in between. 
Slow orbits may be assumed to be static in the time frame of the perturber which leads 
to the impulse approximation. In the opposite extreme, the perturber is nearly static in 
the orbit's time frame which is the adiabatic limit. Most studies appeal to the harmonic 
oscillator model to argue that all orbits with time scales greater than that of the external 
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variation will be invariant; precisely, if the perturbing frequency v is smaller than the 
oscillator frequency then the net change in energy of a random phase ensemble of 
oscillators is exponentially small in their ratio, e.g. proportional to e~^^. Accordingly, 
orbits with < v receive a kick by the perturbation computed using the impulse 
approximation, but the actions of orbits in the adiabatic limit, > z/, are assumed to be 
conserved. The seminal study of globular cluster evolution by Ostriker et al.(1972) was 
based on impulse heating and the adiabatic criterion. This paradigm of a "gravitational 
shock" is now ubiquitous. Chernoff et al.(1986, see also Chernoff and Shapiro 1987) make 
the harmonic model explicit, extending the impulse approximation using Spitzer's (1958) 
treatment of tidal distortions using the linear oscillator model directly. 

This paper presents a theory which allows the evolution to be computed for a general 
time-dependent perturbation regardless of rate. The investigation predicts significant 
evolution in both the adiabatic and impulsive regimes. The explanation for this surprising 
result is discussed in the preceding companion paper (Paper I) which shows that the 
harmonic oscillator model does not apply to the general nonlinear multidimensional orbit 
and describes the failure of adiabatic invariants. To summarize, all stellar systems are 
multidimensional systems with at least two degrees of freedom. Each independent degree 
of freedom has a frequency of oscillation, fij = w r Since any quantity, Q, describing 
a smooth quasi-periodic orbit is well-represented by a rapidly converging Fourier series, 
Q = XjQl ex p(*l • w), the evolution of an orbit is equivalent to the evolution of a set of 
pendula whose frequencies are 1 • fl where 1 is a vector of integers. As long as the frequency 
of the perturbation remains small compared to 1 • f2, each pendulum will be adiabatically 
invariant and the constants of motion (actions) for the original orbit will be conserved. 
However, if one of those frequencies 1 • fl is zero or nearly so, then this term will receive a 
kick from the perturbation. Since one of the Fourier coefficients will change the original 
orbit will no longer conserve its actions. 

For a definite physical example, take a two-dimensional system which has two distinct 
frequencies, fii and fi 2 . If there exists two integers li and / 2 with /ifi — / 2 fi 2 = 0, then for 
every / 2 periods in degree of freedom 1, the orbit will have executed li periods in degree of 
freedom 2, returning to its original configuration. Since the perturbation frequency is small 
compared to either fii and fi 2 , the perturbation "hits" the orbit repetitively at the same 
phase, cumulatively producing a measurable change in the trajectory. This is similar to a 
resonance but with arbitrarily small resonant frequency. Now, a realistic stellar system has 
frequencies continuously distributed in some range. For many orbits, 1 • fl will be far from 
zero for relatively small values of lj and these orbits will be adiabatically invariant. However, 
there will almost always be some combinations of small integers for which 1 • fi = and 
those orbits will not be adiabatically invariant. Averaging over the whole stellar system, 
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the orbits with broken invariants can give an appreciable overall change with magnitude 
similar to the impulsive contribution. 

The large fraction of this paper describes a general method for applying these ideas to 
stellar systems (§2). We begin with a general statement of the time-dependent perturbation 
theory for spherical stellar systems and, for a concrete example, apply the scheme to the 
gravitational shocking of globular clusters by disk passage (§3). Although gravitational 
shocking is the targeted application, this development holds for any finite-duration 
perturbation to a spherical stellar system. The approach could be easily extended to disks. 
The long-term change in energy, action and distribution function are specifically computed. 
The special cases of binary star evolution and very slow perturbations are discussed (§4). 
The effects of the new disk shocking formalism presented here on globular cluster evolution 
is illustrated in Paper III which applies it to Fokker-Planck models. 



2. Problem statement 



Assume a spherical background with the distribution function / (I) where I is the 
action vector. Let the perturbing potential which causes shocking be given by V p (r,t). In 
action-angle variables, the linearized collisionless Boltzmann equation may be written: 

dh dHo_ dh_dh 9H 1 = Q m 
dt 81 ' dw 81 ' dw U 

where Hi = V p and w is the vector of angles conjugate to the action vector, and the 
subscripts '0' and 'I' indicate unperturbed and perturbed quantities and H is Hamilton's 
function. 

This equation may be solved for fi by Fourier transforming equation (1) in angles and 
Laplace transforming in time. The Fourier transform of a phase space quantity Q is given 
by: 

Q = £<2i(iK 1 " w , (2) 

1 

1 ^we-^g^w), (3) 



^ (2tt) 3 

where 1 is a triple vector of integers whose values range over all integers unless otherwise 
noted. Denote the Laplace transform of a quantity Q by Q. Now assume that the 
perturbation vanishes at some time in the past — that is, the cluster used to be far from the 
perturbation — so that f\(t — > — oo) = 0. Then we may solve for f\\ 

i\ ■ dfo/dl 
s + i\- fi(I) 



h = , TZ„ v P i- (4) 
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The perturbed distribution function then follows from the inverse Laplace transform, and 
may be written as 

f 1 = J2Ml,t)e a - w . (5) 
1 

Physically, fi may be thought of as the "wake" in the cluster induced by the perturbation. 

Most of the work in solving for fi is in the step between equation (4) to (5). This 
calculation does not include the self-consistent change implied by the new distribution 
function; this could be done (e.g. Weinberg 1989) but only with considerably more 
computational work. In the next section, we will consider the case of disk shocking and 
perform the inversion explicitly. 



3. Disk shocking 

Let us consider a globular cluster passing through an one-dimensional slab, 
representing the disk. More general perturbations may be done similarly; in particular, the 
time dependent forcing of a cluster on an eccentric orbit will be the subject of a later paper. 

Let v z be the velocity perpendicular to the slab; clearly the transverse velocity is 
irrelevant. To determine the acceleration relative to the cluster center, the perturbing 
potential may be expanded about the cluster center of mass to get the tidal strain: 



1 d 2 K 



V ~ - 

p ~ 2 dz 2 



2 



Z 



(6) 



z c (t) 

where z c (t) is the center of the cluster relative to the slab at time t and z is the vertical 
position in the cluster relative to its center. The truncation error in equation (6) due to 
this tidal expansion is 0([ri/ 2 / h] 2 ). Since ri/ 2 /h ~ 1/60, this error is roughly ~ 0.3%. 

The development below will treat a Gaussian vertical disk profile, 

p(z) = p e-^?lh* (7) 

where h is the scale height of the slab. The cluster center then crosses the disk's midplane 
t = 0. The Laplace transform (two-sided, e.g. van der Pol & Bremmer 1955) for the 
Gaussian profile yields 

2 v z 

The constant quantity V = 4:TrGp h 2 follows from Poisson equation (see Appendix for 
additional details). Other slab profiles may be treated similarly. The exponential slab, 

p(z) = p e-\ z+zM \ /h , (9) 
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will be briefly discussed below for contrast. Although the exponential slab is pathological, 
having a density cusp at z = 0, both models yield similar overall heating rates, remarkably. 
However, the exponential disk does have different asymptotically behavior than an 
everywhere smooth model. 

To complete the calculation defined in §2, we need the expansion of z 2 in action-angle 
variables. Assume the spatial part of V v may be expanded in a harmonic series in spherical 
coordinates: V p (r,6,<f)) = ^2Y nm (6,<f))f nrn (r). Using Tremaine & Weinberg (1984 eqn. 54, 
hereafter TW), each term in the sum may be expanded in action-angle variables: 

oo oo / / 

v P (r,e,<l>) = E E E E Vn^wlU^V 1 ^ ( 10 ) 

/ = li= — OQ l2 = — l ls = — l 

where 

W k nlm = — r d Wl e~^ f nm {r)e U{ ^ (11) 

Z7T J-tt 

and f3 is the inclination of the orbital plane (using the notation of TW). The function 
fi(wi) = if) — w 2 is the difference between the mean azimuthal angle w 2 and the azimuthal 
angle in the orbital plane (cf. TW eq. 38). The fact that m = for all terms above 
demands that / 3 = 0. The quantity V n i m = for |/| > n or |m| > n and non-zero only 
when n + / is even. For example, if n = then / = and the only non- vanishing term is 
VUqqq. Since z 2 = r 2 cos 2 9, f nm (r) = r 2 which is independent of n and m. We may therefore 
simplify the notation for W, defining X 1 ^ = W{] 2 { and write 

X/V 1 -™ (12) 

The inverse Laplace transform of equation (8) may now be straightforwardly performed 
using equations (4) and (12) and retaining only lowest order contributions to get fi(t); 
details are given in the Appendix. 

4. Phase-averaged perturbations 

Most commonly, the long-term it — > oo) change in a moment or the distribution itself 
is desired. The former will be explored in this section. The moment changes are expectation 
values over a specified ensemble. For example, an ensemble of orbits with energy E initially 
will have a spread of energies after the disk passage; the quantity computed in this section 
tells us the average change per orbit in the ensemble. On the other hand, the overall 
evolution is best indicated by the change in the distribution itself. However, equations (4) 
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and (5) show that the phase average of fi vanishes; the second-order contribution, / 2 , does 
not and will be considered in the next section. 

Equation (5) and the explicit form of fi(t) allows us to determine the mean change of 
any constant of the motion using Hamilton's equations. For any function of phase space, 
we have (e.g. Goldstein 1950) 

dQ__dQ_ 
dt ~ dt 

where [] are Poisson brackets. The first term on the right-hand-side vanishes if Q is a 
constant of the motion. On expanding the previous expression, the lowest order terms 
vanish by assumption. Retaining only first-order terms yields 

dQ dQ dH! d#i dQ 



+ [Q,H] 



(13) 



dt dw dl 

The angle or phase average of Q is then: 

1 



dw dl 



(2,Y J dw d ft h = £ 



dQ dV pl 
dw dl 



dQ 
dl ' pl 



In particular, for a component of action we have 



iAV) = E-Wi(M) = J2^v P -iMi,t) 
i i 



(14) 



(15) 



(16) 



and since E = H n (l) we have 



i i 



(17) 



Finally, the phase-averaged change in Q for an ensemble of fixed action I may then be found 
by integration: (AQ) = dt (Q(I,t)). 



4.1. Change in energy 



As an example, let us evaluate (AE) explicitly. For the spherical system discussed 
here, Ii is the standard radial action, I 2 = J and I 3 = J z where J and J z are the total 
and z projection of the angular momentum. If the system has no net rotation then the 
distribution function will be independent of J z and we may integrate over I 3 keeping I 2 
fixed. This motivates changing variables from I 3 to inclination f3 defined by cos f3 = l 3 j I 2 - 
Then, dl = dl^d^^dfi sin f3 and we may define: 

((Q)) = f J/3sin/3(Q). (18) 
Jo 
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We now evaluate ((AE)) resulting from the passage through the slab using the 
expression for f\ and V p \{t) and find (see Appendix) 

\\ 1 1 4 ft 4 V V15 90 2 J 1 ' 2 ' 



x l-#44L^e- fc2(1 ' n)a/ ^. (19) 
<9I 2(v z /h) 2 K ' 

The analogous expression for the exponential disk is also derived in the Appendix: 
KAF)) - ( 2l ) 3y »^ f 1 + 43 A WP 
j a/ 2(l-n)K/tf 

91 [(^//,)2 + (l.fi)2]2- ^ J 

If the distribution function depends on energy alone, f = f(E), then 

1^-1 <9#4/o_, / T ^ / 91 N 

1 ' di " 1 ' di dE " ' ai " ' ( W ( j 

Upon substituting into the expression for ((AE)) } we find that to lowest order, disk 
shocking increases the energy for all E if df /dE < 0. 



4.2. Comparison of Gaussian and exponential disks 



Features of the expressions for ((AE)) are worth noting. Both equations (19) 
and (20) scale explicitly as ((AE)) oc [(1 • Sl)/(v z /h)] 2 in the impulsive limit 
[(1 • SI) / (v z / h) <C 1] as expected. However in the adiabatic limit, equation (20) is 
proportional to [(1 • SI) / (v z / h)]~ 2 while equation (19) appears to be exponentially small. 
However, if there are commensurabilities between frequencies, 1 • SI = 0, the power law 
behavior does obtain for equation (19) as well for the physical reasons described in Paper I. 

To see this, take an ensemble at fixed energy but isotropic in velocity. Integrating the 
distribution over the ensemble requires an integration over dJJ/fti(E } J). The integrand in 
this case is proportional to 

(1 • SI) 2 h 2(\.Sl ) 2l2vl 

2(v z /h) 2 

The variation of the exponent depends most strongly on the variation of the ratio of 
frequencies \ = ^2/^1 for small v z /h. In these terms, the integral to be done is roughly 
/ ^XX 3e_c * 2 ' x ~ x °' 2 where the weaker dependence of (A^ 1 ! 2 on \ nas been ignored. If the 
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variation of \ includes \o, the integral gives a contribution proportional to a 4 and leads to 
the scaling ((AE)) oc [(1 • n)/(v z /h)]~ 2 in the adiabatic limit [(1 • £l)/(v z /h) > 1]. Explicit 
numerical evaluation shows that the Gaussian case does have a power law dependence but 
with steeper slope. Figure 1 compares ((AE)) for the Gaussian and exponential slabs. Both 
models yield very similar results suggesting that the response to the gravitational shock is 
weakly dependent on the disk's vertical structure 

Conversely, the reader may be surprised that equation (20) shows no exponential cutoff 
for incommensurate orbits in the adiabatic regime. This is a consequence of the exponential 
disk profile, whose force is not smooth at z = 0. This discontinuity provides a power at all 
temporal frequencies, and loosely speaking, creates resonant interactions at all frequencies. 

4.3. Check by simulation 

Figures 2-4 compare the theoretical predictions of equation (f 9) with restricted n-body 
simulations for a Wo = 5 King model. The integration accuracy corresponds to an absolute 
error in energy of approximately 10~ 7 per orbit. The parameters roughly describe a typical 
globular and Galactic disk — M = 3 x fO 5 M , r t = fOOpc, R g = 8kpc, V c = 200 kms" 1 , 
and h = 325 pc. However in order to improve signal/noise, the disk potential is artificially 
increased by a factor of 3 over its realistic value; this enables the resolution of the highest 
frequency encounter (Fig. 4). The bin with largest energy is undersampled due to escaping 
stars and the bins with the smallest energy are affected by the integration accuracy, 
indicated by the dashed horizontal line. The agreement is good for energies between these 
two limits. 

4.4. Application to binary stars 

The gravitational shocking theory may be used to compute energy changes for an 
ensemble of binary stars and is consistent with the standard results obtained for binary 
star systems (Heggie 1975). In two-body problem, the orbital frequencies are degenerate, 
fii = fi 2 an d ^3 = 0, which means that 1 • fl vanishes if and only if h vanishes. Therefore, 
the contribution to ((AE)) will always vanish if 1 • fi = and the contribution will be 
exponentially down in the ratio (1 • Sl)/(v z /h) as expected for 1 • fi ^ 0. The energy is a true 
adiabatic invariant for binary stars. 

However, the behavior is different for the changes in individual actions ((AI 3 )). 
Equation (16) gives the change in a single action, Ij and ((AI 3 )) may be derived in 
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the same way as ((AE)). The distribution function for a single orbit with action I is 
/(I) = 8(1 — I )/(27r) 3 . Integrating by parts and using the Gaussian disk perturbation for 
definiteness yields: 

{{AIi)) = iV^JW^ Sh0 ll5 + 90 H lil ' di 

In the limit Cl/(v z /h) ^> 1, terms with 2 (? 7^ will be exponentially down, but unlike the 
case for AE } the terms with ^Ij = do not vanish. In fact, in this limit the previous 
equation becomes 

Evaluating (A^p explicitly shows that this expression for ((AI 3 )) does not vanish (except 
for purely radial or circular orbits) and drives orbits to larger eccentricities even though 
energy is adiabatically conserved (((A7i)) = ((A7 r )) > and (((A7 2 )) = ((AJ)) < 0). It is 
conceivable that this mechanism modifies the primordial elements of young binaries in the 
dense environments in which they form and in longer lived star clusters. 




4.5. Evolution for very slowly changing perturbations 



Let us take a general perturbation V^(x, t) which is turned on in the past and turned off 
in the future. Now if the time variation is made very slow, e.g. V v = V^(x,t/r) for large r, 
the scaling for the long-term change in the distribution may be computed directly (originally 
suggested and computed by Scott Tremaine 1994). The calculation is straightforward and 
is sketched for ((AE)) in the Appendix. Integrating over the phase space for the entire 
model gives the total energy change for the system due to a gravitational shock: 



«<A£))) 



I 



4- 3 EE 



dss 2 du n 



I*; 



. (24) 

where ^i\ 2 ; (I, s) = Vi \ 2 \ z (f3)W l j 2 ls (I) is the action-angle transform of the perturbation 
which has been Laplace transformed in time and 8 is the Dirac delta function. Equation 
(24) shows that the overall change in energy for a cluster or dwarf galaxy scales as I/r. 
Comparison with equation (19) and the derivation in the Appendix suggests that I/r 
scaling from the peak amplitude at r I/O, should be roughly correct. At any rate, this 
is very weak adiabatic cutoff, qualitatively and quantitatively different than the expected 
invariance at large r. 
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This scaling shows that no realistic stellar systems are truly invariant to slowly 
changing external perturbations. For example, the core of a galactic cluster will be modified 
by the cluster's vertical oscillations even though the oscillation period may be many times 
the stellar orbital periods. In addition, it is often assumed that a cannibalized dwarf with 
phase space density higher than its captor will survive disruption. However, this result 
suggests that gravitational shocking may disrupt even dense dwarf companions. 



5. Second-order calculation of perturbed distribution function 



Here, we will compute the overall change to the distribution function directly which 
requires extending the calculation to second order. The perturbed distribution function 
allows the long-term evolution of the galaxy or star cluster to be estimated. This result will 
be used in the Fokker-Planck implementation in Paper III. The features of the solution are 
explored below using the singular isothermal sphere. 



5.1. Derivation 



The first-order distribution function was solved in §2. The Boltzmann equation to next 
order equation is: 

d^ + d^ m d^ + dH lm df L _df Lm dH 1 = Q 
dt dl dw dl dw dl dw 

The H 2 term vanishes in the absence of self-gravity, Following §2, the Fourier-Laplace 
transform of the second order equation becomes 

— -e ■2^i\fi„e -2^^T e -l^iWye \ = 

J; ,JS - in in ,JS - J; ) 

(26) 

where f 2 are the expansion coefficients of the second-order distribution function. We 
may restrict our attention to the 1 = term which gives the only contribution to the 
phase-averaged secular change. The integral in w may then be done immediately, yielding 
the solution 

1 poo f) 

/ 2 l=o = — / dte- st Y.iV -^{V-vfv). (27) 

S J -co j, Ol. 
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Using the action-angle expansion of the Gaussian disk perturbation, the time integral may 
be performed explicitly. 

We must now do the inverse Laplace transform to get from i*2l=o(I? s ) to i*2l=o(I? t)- The 
joint domain of convergence is > and we deform the contour to the imaginary axis 

to perform the inverse transform for t — > oo. Since there is a pole along the imaginary axis, 
we may divide the contribution in to the principal part and residue. Changing variables 
to s = iy shows that the principal part does indeed exist. However, it has a factor of the 
form lim^oo exp(iyt). Because the rest of the integrand is analytic on physical grounds, the 
principle part vanishes. The final expression consisting of the residue contribution alone is: 

fwl =°- ^2(v z /hy^J di\ e 1 di [zh[zU j {Z * } 

having used the symmetry in 1 • over the summation in 1. The quantity [z 2 }\ is the 1 
Fourier component of the action-angle expansion of z 2 whose expansion is given in §3. 
Integrating over all orbital planes this becomes: 

it \ - V Q % V* x f 1 43 k \^ d S 
{fwl=o) ~ 4^2(^7/0V^ Hl5 + 90 '^'all 



-(l.fi) 2 /2(^// l ) 2 , . ^Jo 
31 



dfo vU 2 | (2g) 



x) 1 



Finally, we may write the distribution function in terms of E and J since this is the most 
common form with 

With these conventions, equation (29) becomes: 

It \ V o * ^ c f 1 7 c \ A o 9 1 9 \ \ 

e -W/W)»|i.n| + i^j^|j. (30) 

In practical applications, (f 2 ) must be evaluated numerically. Most of the time in evaluating 
(f 2 ) and ((AE)) is the computation of (A^l 2 . An easily implemented numerical approach 
is described in the Appendix. 
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5.2. Isothermal model 



The singular isothermal sphere is similar to the King model over part of its energy 
range but because it is scale free, the disk-shocking calculation need only be done for a 
single energy over the range of desired v = v z /h. The simplicity of this model makes it a 
good tool for exploring the basic features of the theory. 

The singular isothermal sphere may be defined as follows. If one takes the potential to 

be 

U(r) = 2<j 2 lnr (31) 

then from Poisson's equation 

M(r) = ^r, (33) 



and the distribution function is 



normalized consistently with p(r). However, we will take f(E) = exp(—E/cr 2 ) for 
convenience in the calculations below. 

Using these relations, we may define scale-free quantities. Defining n = J / J max (E) } 
where J m ax(E) is the maximum orbital angular momentum at fixed energy E, it is easy to 
show that 

Jraax^E) Jmax(^)& ^ •> 

r(E,K) = r(0 } K)e E ^\ 

and 

nj(E,K) = n J (o,K)e- E/2a \ 

The quantity r is the radius of an orbit with E and n at a particular phase, e.g. a turning 
point. Substituting these expressions into equation (30) shows that (f 2 ) is independent of 
energy except through the factor exp[— (1 • £l) 2 /2v 2 ]. Therefore, (f 2 ) as a function of E for 
v = 1 is equivalent to (f 2 ) evaluated at E = for v = exp(£ , /2cr 2 ). This is also true for 
((AE)). Other quantities may be scaled similarly. 

Figure 5 shows the change in energy in an ensemble where orbits have single value of 
E initially, ((AE)). Note that this is not the same as the change in energy of the ensemble 
with post-shock energy E. The top panel is multiplied by v 2 to show the asymptotic 
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behavior in the impulsive regime, ((AE)) oc v~ 2 . The lower panel is scaled to show ((AE)) 
at E = as a function of lnz/ 2 . Detailed analyses show that ((AE)) oc z/ 3 ' in the adiabatic 
regime (E <J — 1). 

The lower panel also indicates significant change in energy per orbit inside of the 
adiabatic boundary. The features in this figure are due to individual contributions (l\ } l 2 ). 
The individual contributions are shown separately in Figure 6. The terms (1,0) and (0,2) 
contribute the peak on the right hand side and comprise most of the impulsive contribution. 
These components have incommensurate frequencies and similar profiles with the lowest 
order having the highest amplitude. The peak on the left hand side is dominated by the 
(1, —2) term and is largely in the adiabatic regime; this contribution is due to the accidental 
degeneracy as discussed in Paper I. The (2, —2) term is also in this regime but is of lesser 
importance. Note the different profiles for each of these and the incommensurate terms. 

Figure 7 shows v 2 (f 2 ) as a function of E / 'a 2 for v = 1 and (f 2 ) at E = as a function 
lnz/ 2 . The value of is approximately 2.5 for E = which corresponds to E = In v 2 2 at 
the boundary between the adiabatic and impulsive regimes, fi/z/ 1. Again, the asymptotic 
impulsive behavior (f 2 ) oc v~ 2 is obtained for E ^ 2. The overall peak contribution is in 
the transition between the adiabatic and impulsive regimes. Contributions of individual 
terms (h 7 l 2 ) dominating (f 2 ) is shown in Figure 8. The incommensurate terms (0,2) and 
(1,0) contribute in the impulsive regime; they are effectively zero for E <^ and decrease 
quickly inside of E ~ 2. The commensurate (1, —2) and (2, —2) terms contribute at smaller 
binding energies. The (1, —2) term is the strongest of these and contributes positive density. 
Responses with positive pattern speed [all but (1, —2)] shift stars outward to higher energy 
and vice versa. Clearly the relative amplitudes of individual contributions depends on the 
particular distribution and therefore the details of the isothermal model are suggestive only. 
As a function of passage frequency, the lower panel in Figure 7 shows that (f 2 ) increases 
(decreases) in the adiabatic (impulsive) regime. 

Since the singular isothermal sphere has infinite extent, it is not meaningful to compute 
the total heating (because of the semi-infinite interval in energy space in the impulsive 
regime). For a truncated cluster, the net effect of the shock on the cluster depends on 
the truncation energy relative to the peak in energy transfer. We may use the singular 
isothermal sphere as a guide to the possible effects of a gravitational shock. First, if peak 
feature in Figures 5 and 7 occurs at relatively large binding energies, the overall contribution 
will be dominated by the impulsive interactions. However, if the feature occurs in the halo 
(e.g. small z/), heating may be dominated by adiabatic interactions. Secondly, although the 
interaction increases the energy of an initially monoenergetic ensemble (cf. Fig. 5), the 
r.m.s. change is large compared to the net change. In other words, individual particles are 
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both gain and lose energy. If the peak occurs at sufficiently low energy in the halo, the 
heated orbits may escape leaving a cluster with overall higher binding energy. It is therefore 
possible to "cool" the cluster in a disk shock. 

6. Thick-disk and bulge shocking 

Thick-disk shocking is a straightforward variant and somewhat easier to compute than 
the calculation above. Since the perturbation has the form V v = g(t)z 2 } if the cluster is 
dynamically part of the disk then git) will be a periodic function of time and may be 
expanded as a Fourier series in its vertical oscillation period, P: 

oo 

g(t) = £ 9ke^\ (35) 

k= — oo 

where u = 2tt/P. The Laplace transform of a Fourier series is trivial and in most cases, 
will converge rapidly with increasing \k\. 

The calculation is analogous for bulge shocking but the potential perturbation 
expansion will be more general than the z 2 dependence and include all second order 
moments (Y 2m terms). In addition the Fourier expansion of git) will have two indices 
corresponding to the radial and azimuthal motions of the cluster orbit. This approach can 
easily include the centrifugal potential but not the velocity-dependent Coriolis force. 

External heating will most likely play a different role in evolving clusters which 
are kinematically a halo component or a disk or thick-disk component. Eccentricity 
is likely to critical also, as Aguilar et al.(1988) have pointed out for bulge shocking. 
A detailed investigation of gravitational shocking for a globular cluster on a general 
orbit — implementing disk, orbit and bulge shocking — is in progress. 

7. Summary 

Paper I showed that adiabatic invariants are NOT exponentially controlled 3 for all 
orbits if the number degrees of freedom for the system is greater than one. Consequently 
for a stellar system, some orbits will be strongly perturbed even if the characteristic 
frequencies are much larger than the perturbing frequency v. This leads to measurable 



3 proportional to e u l v where u> and v are the characteristic system perturbation frequencies 
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overall heating ol a star cluster or galaxy in the adiabatic regime. The change in energy of 
an initially monoenergetic ensemble scales approximately as (u/fl)^ (2^/3^3 depending 
on the model) compared to (z//fi)~ 2 for the impulsive regime; the two connect smoothly at 
v 0. For a perturbation which changes slowly over time r, the total change in kinetic 
energy due to the gravitational shock scales as 1/r; there is no sharp adiabatic cutoff. 

Heating in the adiabatic regime requires frequencies that are irrational ratios of each 
other except on discrete surfaces in phase space. 4 This is a very weak condition, true for 
nearly all commonly used cluster and galaxy models. A binary star system is counter 
example since fii = fi 2 7^ 0, fi 3 = for all energies. Nevertheless, even though the energy 
of a binary is adiabatically invariant, the method shows that individual actions for a 
binary system are not invariants. For example, a slow perturbation can change a binary's 
eccentricity in the adiabatic regime. 

In addition to energy changes, we derived the long-term change to the distribution 
itself which requires extending the perturbation theory to second order. The expression for 
the shocked distribution is easily incorporated into a Fokker-Planck simulation, and this is 
the subject of Paper III. 

Since the theory predicts gravitational shocking over a much wider range of encounter 
rates and at larger magnitude than previous estimates, a variety of scenarios in addition 
to standard disk shocking may require revision. First, the adiabatic criterion does not 
abruptly limit the heating of clusters due to the relatively slow vertical motion in the disk. 
The development may be modified to account for a periodic perturbation appropriate to 
shocking a thick-disk globular cluster and by cumulative encounters with GMCs. Secondly, 
the time-dependent external force of a globular on an eccentric orbit is also a gravitational 
shock, often called a "bulge shock" in its extreme form. This has already been discussed 
by Aguilar et al.(1988) among others but the current work suggests that its importance 
may be even greater. The same effects may be important for dwarf galaxies on moderately 
eccentric orbits and may be easily applied to disks. This approach will not account for 
the Coriolis force but will be correct for a radial cluster orbit and detailed simulations 
(Murali et al.1994) suggest the Coriolis force causes no dominant effects. The seminal 
shocking problem, the response of a star cluster to a passing molecular clouds, a point-mass 
fly by (Spitzer 1958) may be performed similarly to include the additional heating. The 
development for binary star evolution is sketched in §4.4 and may be useful for computing 
the effects of protostellar evolution. 



4 commensurabilities 1 • ft = where 1 a vector of integers 
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A. Gaussian slab profile 



The potential for the Gaussian slab, p(z) = p exp(—z/h), is 



(m; 



where V = 4:TrGp h 2 . Using the tidal prescription in §3 with z = z c (t) = v z t } the (two-sided) 
Laplace transform is 



2 v z 



(A2) 



The inverse transform required to derive the distribution function (cf. eq. 5) may be 
performed by deforming the contour to the imaginary axis, since the joint domain of 
convergence is 3?(s) > 0. Simplifying, the required integral takes the form 



Z(t) 



1 /*°° 



dy 



1 



(A3) 



2tt J-oo j i(y + 1 • 

We now separate the integral to the residue and principle part denoted by TZ and V 
respectively. The residue trivial. The principle part may be computed by defining 
z = y + 1 • f2, splitting the integrand into the two semi-infinite intervals ( — oo, 0] and [0, oo), 
changing z — > —z in the former term and combining. Putting these together with equation 
(4) yields 

1 r°° dz 



-h 2 z 2 /4v 2 z 



TTl JO Z 



smh — „ , , „ cos tz — % cosh — „ , , „ sin tz 



2v 2 Jh 2 



2v 2 Jh 2 



(A4) 



The expression [z 2 }\ denotes the 1 Fourier coefficient in the action-angle expansion of z 2 (eq. 
12). 

We now evaluate ((AE)) for the using the expression for f\ and V p \{t). In particular, 
we need to do the integral 



Y 



dte h 2 (t+ Zo /v z ) 2 /vl^ n + V y 



(A5) 
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Performing the algebra, one finds 

-h 2 (l-Sl) 2 /2v 2 z 



Y 



1 + 



dz 



-W* sinh 



TTl JO Z 



i • n 



(A6) 



From its definition in TW, we have 

V n i 2 iM) = ^ 2 / 3 (/3)K/ 2 (vr/2,0)z' 3 -' 2 . Using the 
orthogonality of the rotation matrices, f£ d/3 sin f3rf 2 ^ (/3)r^ 3 = 2/(2/ + l)6 nn i, and the 
relations K-/-m = ( — l) m K/m and W k a * m = ( — l) m W~~_^_ m we see that only terms with the 
same value n will contribute. As in §3, only the terms even in 1 • fi contribute to ((AE)) } 
which eliminates the second term in Y. So finally, the expression becomes (eq. 19): 



«A£)) 



h 4 



1 



43 



■TZ>3 (^ + ^ 2 oj |A /2 



. a/ i-n _ 



15 90 

h 2 (l-Sl) 2 /2v 2 z 



h |2 



(A7) 



<9I 2(^//i) 2 

To minimize the computational expense, the sum over 1 may be written as 

E = 2 £' • 

1 '1=0 

i 2 = -2,0,2 
i 3 = -oo 

The restriction in the limits of the sum follow from the properties of Y n \ m and the prime is 
to remind that l 2 ^ — 2 if li = 0. 



B. Exponential slab profile 



The potential of the disk profile, 

p(z) = p e 

follows from direct integration of Poisson's equation yielding 



-\z+v z t\/h 



V(z) = 47rp h 2 



h 



+ e~^ h - 1 



:bs) 



Note that the third derivative of the potential or second of the force is discontinuous (see 
text of a discussion of the implications). 



The Laplace transform is 



z 2 V 
V = — — 

9 2 V 



1 1 

+ 



v z /h — s v z /h + s 



(B9) 
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with the domain of converge —v z /h < < v z /h. Using equation (4), we now take the 

inverse Laplace transform to get fi(t). Since the domain of convergence for f\ is > 0, 

the joint domain is < < v z /h and since the poles are explicit, the inverse is simple. 



To facilitate evaluation, the contour may be deformed to c 
t > t„. One finds 



oo if t < t„ and c 



-oo 



if 



fl(t < tc 



fl(t > tc 



df^Vo_ 
dl 2V 

dhVo_ 
dl 2V 



v z t/h 



v z /h + i\ ■ f2 



2v z /ht 



il-Slt 



-v z t/h 



A 2 



,/h 



i\ ■ n 



;bio) 



where A 2 = (v z /h) 2 + (1 • fi) 2 . 



Proceeding as in the previous section, the time integration may be done by breaking the 
interval into the two parts ( — oo, 0] and [0, oo). Again since only terms even in 1 contribute 
by symmetry, only one term contributes. Integrating over f3 then yields (eq. 20) 



«A£)) 



(27T? V 



2 oo 



h 



1 • 



df 2(\-n)(v z /h) 
ai [(v z /hy + (i • o) 



/ 1 43 

V15 ' 90 5 ' 2 ' 

2 



lE^30(^ + ^ 2 0)|^| 2 



t2 - 



:B121 



C. ((AE)) for very slow perturbations 



In order to show ((AE)) is proportional to 1/r as r — > oo, the same steps in the 
previous two appendices are repeated formally, expanding at the end in powers of t/r. In 
particular to start, we need to do the integral 

/oo 
V p -i(t)Mt), (C13) 
-oo 

where fi(t) is given by the inverse transform of equation (4) and the perturbed potential 
is V p (x.,t) = U(x.,t/r). Assuming that the Laplace transform for V v converges at least for 
Re(s) > 0, substituting the explicit expressions for inverse transforms into equation (C13) 
and using Cauchy's theorem gives 

Z = [°° dtU(x,t/T) [°° dt'e-^-^Q^-t'p^t' '/t) 

•J -co J — CO 

= r dqe- thnq r dtU(x,t/T)U(x,(t + q)/T) (C14) 

JO J-oo 
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where denotes the Heavyside function and the second equality follows from the change of 
variables q = t — t' followed by t — > t + q. Finally, using the convolution theorem and the 
even symmetry of the result in 1 gives 



Z = -\U(^\.SlT)\\ 



(C15) 



Now consider the expression 



dL0L0 2 \U{^L0T)\ 2 f{t0) 



for some arbitrary function f{oo). In the limit r — > oo, the quantity |[/(x,cu)| 2 represents a 
positive definite distribution in u peaked at u = 0; this may be represented as a Fourier 
transform in u. Using this explicitly and after some manipulation, one finds 



duu 2 \U{^u)\ 2 f{u) = -f(0) I dss 2 \U(^s)\ 



and therefore in the large r limit 



lo 2 \U{^lo)\ 2 = -6{lo) / dss 2 \U(^s)\ 



Similarly, one also finds 

L0\U(^L0)\ 2 

Putting this together gives 



S'(uj) / dss 2 \U(^s)\' 



Z T j 



i • n^s(\ ■ o) - ij^8'{\ ■ n) 

dE y ' d.J y ' 



(C16) 



(C17) 



(C18) 



dss 2 \U l {l,s)\ 2 . (C19) 



D. Computational considerations 



Calculation of the central equations of this paper, equation (30) in particular, should be 
straightforward to anyone familiar with solving for cluster evolution using the Fokker-Planck 
approach. The only new detail is the computation of the potential transform X 1 ^ . The 
following procedure solves for all needed quantities by direct quadrature and a set of 
coupled ordinary differential equations. 

1. Given E and J, determine the turning points r a and r p by solving the standard 
implicit equation = 2[E — U(r)] — J 2 /r 2 . 
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Table 1: System of equations for computing potential transforms 



Variable 



Equation 



Initial condition 



V = 
f 

Wi 

x! 1 

<2 



y 



y = 
f = 

Wi 



dU J 2 



X 



Ox 
Ox 

7T 



cos(/iu; — hf)^ 2 



r(t = 0) = r p 
y(t = 0) = 

f(t = 0) = 
Wl (t = 0) = 

xNt = o) = o 



2. Determine the orbital frequencies by direct quadrature: 

1 



7T 


/ dr 






o 2 = 








7T 7r p 



^2[,E - £/(r)] - J 2 /r 2 ' 

1 J 



^/2[E - U(r)] - J 2 /r 2 r ^ 

(D20) 



3. The potential transform may now be done by simultaneously integrating the radial 
component equations of motion for the orbit, the differential expressions for the radial 
angle, the difference between the true and mean azimuthal angles, and the potential 
transform itself. The set is shown in Table 1. The integration is assumed to begin 
at apocenter at t = 0. For a particular pair (7i,/ 2 ), the system has five coupled 
equations. However, since only a small number of terms /i, / 2 contribute as shown in 
the previous section and §5, all n of these may be done simultaneously, giving a set of 
4 + n equations. 
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FIGURE CAPTIONS 

Fig. 1. — Comparison of ((AE)) for a Wo = 5 King model perturbed by a Gaussian slab 
(solid) and exponential slab (dashed) for v z /h = 0.03,1.0,10.0 (labeled at right). The 
amplitude is arbitrary. 

Fig. 2. — Comparison of direct simulation to perturbation theory calculation. The histogram 
shows the results of ((AE)) computed from a direct integration in a fixed potential with 200K 
particles realized from a Wo = 5 King model. The dotted line shows the mean integration 
error per orbit in each bin. The solid curve shows the predicted relation using the formula 
in the notes. 

Fig. 3. — Same as Fig. 2 but with the passage frequency decreased by a factor of 3. 

Fig. 4. — Same as Fig. 2 but with the passage frequency increased by a factor of 10. 

Fig. 5. — The energy change for ensembles whose orbits have the same value E initially. 
The top panel shows ((AE)) multiplied by v 2 but for true E and the bottom panel scaled 
to E = 0. 

Fig. 6. — As in the lower panel of Fig. 7 but but individual contributions (/i, / 2 ) are shown. 
The total is omitted for clarity. 

Fig. 7. — The second-order perturbed distribution function multiplied by v 2 = (v z /h) 2 (top 
panel) and scaled to E = as a function of In v 2 (lower panel). The lower scale may used to 
read v 2 (f2) as a function of £ or as a function of v at E = 0. 

Fig. 8. — As in the upper panel of Fig. 7 but individual contributions (h^h) are shown 
along with the total. Other terms are negligible. 



